In [ ]:
from __future__ import division, print_function
import os
import sys
from six.moves import cPickle as pickle
# Third-party
import astropy.coordinates as coord
import astropy.units as u
import matplotlib as mpl
import matplotlib.pyplot as pl
%matplotlib inline
pl.style.use('apw-notebook')
import numpy as np
# Custom
import gala.dynamics as gd
import gala.integrate as gi
import gala.potential as gp
from gala.units import galactic
import ophiuchus.potential as op
from ophiuchus.data import OphiuchusData
from ophiuchus.util import integrate_forward_backward
from ophiuchus.coordinates import Ophiuchus
from ophiuchus import galactocentric_frame, vcirc, vlsr, RESULTSPATH
from ophiuchus.experiments import MockStreamGrid
from ophiuchus.util import get_potential_stream_prog
from ophiuchus.plot import surface_density
plotpath = "/Users/adrian/projects/ophiuchus-paper/figures/"
if not os.path.exists(plotpath):
os.mkdir(plotpath)
In [ ]:
ophdata = OphiuchusData()
ophdata_fit = OphiuchusData("(source == b'Sesar2015a') | (Name == b'cand9') | (Name == b'cand14')")
ophdata_fan = OphiuchusData("(source == b'Sesar2015b') & (Name != b'cand9') & (Name != b'cand14')")
In [ ]:
all_names = ["static_mw"] + ["barred_mw_{}".format(i) for i in range(1,10)]
short_names = ["static"] + ["bar{}".format(i) for i in range(1,10)]
name_map = dict(zip(all_names, short_names))
In [ ]:
# FOR PAPER
style = {
"linestyle": "none",
"marker": 'o',
"markersize": 2,
"alpha": 0.1,
"color": "#555555"
}
line_style = {
"marker": None,
"linestyle": '--',
"color": 'k',
"linewidth": 1.5,
"dashes": (3,2)
}
data_style = dict(marker='o', ms=4, ls='none', ecolor='#333333', alpha=0.75, color='k')
data_b_style = data_style.copy()
data_b_style['marker'] = 's'
In [ ]:
# contour stuff
levels = np.array([-2,-1,0,1]) - 0.1
for k,name_group in enumerate([all_names[:5],all_names[5:]]):
fig,allaxes = pl.subplots(3, 5, figsize=(9,6.5), sharex='col', sharey='row')
for i,name in enumerate(name_group):
axes = allaxes[:,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot, streams, prog, _ = get_potential_stream_prog(name)
stream = streams[:,best_ix]
model_c,model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# cut to just points in the window
ix = (model_c.l > 0*u.deg) & (model_c.l < 11*u.deg) & (model_c.b > 25.5*u.deg) & (model_c.b < 34.5*u.deg)
# density in sky coords
sky_ix = ix & (model_c.distance > 4.5*u.kpc) & (model_c.distance < 10*u.kpc)
grid,log_dens = surface_density(model_c[sky_ix], bandwidth=0.2)
cs = axes[0].contourf(grid[:,0].reshape(log_dens.shape), grid[:,1].reshape(log_dens.shape), log_dens,
levels=levels, cmap='magma_r')
axes[1].plot(model_c.l[ix], model_c.distance[ix], rasterized=True, **style)
axes[2].plot(model_c.l[ix], galactic.decompose(model_v[2][ix]), rasterized=True, **style)
# l,b
axes[0].plot([3.8,3.8], [31.5,32.5], **line_style)
axes[0].plot([5.85,5.85], [30.,31], **line_style)
axes[0].set_aspect('equal')
# l,d
axes[1].plot([3.8,3.8], [8.7,9.3], **line_style)
axes[1].plot([5.85,5.85], [7.2,7.8], **line_style)
# l,vr
axes[2].plot([3.8,3.8], [276,292], **line_style)
axes[2].plot([5.85,5.85], [284,300], **line_style)
# the data
# _tmp = data_style.copy(); _tmp.pop('ecolor')
# axes[0].plot(ophdata_fit.coord.l.degree, ophdata_fit.coord.b.degree, **_tmp)
_tmp = data_b_style.copy(); _tmp.pop('ecolor')
axes[0].plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **_tmp)
# axes[1].errorbar(ophdata_fit.coord.l.degree, ophdata_fit.coord.distance.to(u.kpc).value,
# ophdata_fit.coord_err['distance'].to(u.kpc).value, **data_style)
axes[1].errorbar(ophdata_fan.coord.l.degree, ophdata_fan.coord.distance.to(u.kpc).value,
ophdata_fan.coord_err['distance'].to(u.kpc).value, **data_b_style)
# axes[2].errorbar(ophdata_fit.coord.l.degree, ophdata_fit.veloc['vr'].to(u.km/u.s).value,
# ophdata_fit.veloc_err['vr'].to(u.km/u.s).value, **data_style)
axes[2].errorbar(ophdata_fan.coord.l.degree, ophdata_fan.veloc['vr'].to(u.km/u.s).value,
ophdata_fan.veloc_err['vr'].to(u.km/u.s).value, **data_b_style)
# Axis limits
axes[0].set_xlim(9.,2)
axes[0].set_ylim(26.5, 33.5)
axes[1].set_ylim(5.3, 9.7)
axes[2].set_ylim(225, 325)
# Text
axes[0].set_title(name_map[name], fontsize=20)
axes[2].set_xlabel("$l$ [deg]", fontsize=18)
if i == 0:
axes[0].set_ylabel("$b$ [deg]", fontsize=18)
axes[1].set_ylabel(r"$d_\odot$ [kpc]", fontsize=18)
axes[2].set_ylabel(r"$v_r$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
fig.tight_layout()
# fig.savefig(os.path.join(plotpath, "mockstream{}.pdf".format(k)), rasterized=True, dpi=400)
# fig.savefig(os.path.join(plotpath, "mockstream{}.png".format(k)), dpi=400)
In [ ]:
import matplotlib.colors as colors
def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
new_cmap = colors.LinearSegmentedColormap.from_list(
'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
cmap(np.linspace(minval, maxval, n))[::-1])
return new_cmap
x = np.linspace(0, 1, 128)
pl.scatter(x,x,c=x,
cmap=truncate_colormap(pl.get_cmap('magma'), 0., .9))
custom_cmap = truncate_colormap(pl.get_cmap('magma'), 0., 0.9)
In [ ]:
scatter_kw = dict(cmap=custom_cmap, alpha=0.8, s=4, vmin=-1, vmax=0, rasterized=True)
data_b_style = dict(marker='s', ms=4, ls='none', ecolor='#31a354', alpha=1, color='#31a354')
line_style = {
"marker": None,
"linestyle": '-',
"color": '#888888',
"linewidth": 1.5
}
In [ ]:
for k,name_group in enumerate([all_names[:5],all_names[5:]]):
fig,allaxes = pl.subplots(3, 5, figsize=(10,6.5), sharex='col', sharey='row')
for i,name in enumerate(name_group):
axes = allaxes[:,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot, streams, prog, release_t = get_potential_stream_prog(name)
release_t = release_t/1000.
stream = streams[:,best_ix]
model_c,_model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# remove progenitor orbit
model_c = model_c[1:]
model_v = []
for v in _model_v:
model_v.append(v[1:])
# cut to just points in the window
ix = (model_c.l > 0*u.deg) & (model_c.l < 11*u.deg) & (model_c.b > 25.5*u.deg) & (model_c.b < 34.5*u.deg)
# plot model points
cs = axes[0].scatter(model_c.l[ix], model_c.b[ix], c=release_t[ix], **scatter_kw)
axes[1].scatter(model_c.l[ix], model_c.distance[ix], c=release_t[ix], **scatter_kw)
axes[2].scatter(model_c.l[ix], galactic.decompose(model_v[2][ix]), c=release_t[ix], **scatter_kw)
# l,b
axes[0].plot([3.8,3.8], [31.5,32.5], **line_style)
axes[0].plot([5.85,5.85], [30.,31], **line_style)
axes[0].set_aspect('equal')
# l,d
axes[1].plot([3.8,3.8], [8.7,9.3], **line_style)
axes[1].plot([5.85,5.85], [7.2,7.8], **line_style)
# l,vr
axes[2].plot([3.8,3.8], [276,292], **line_style)
axes[2].plot([5.85,5.85], [284,300], **line_style)
# the data
# _tmp = data_style.copy(); _tmp.pop('ecolor')
# axes[0].plot(ophdata_fit.coord.l.degree, ophdata_fit.coord.b.degree, **_tmp)
_tmp = data_b_style.copy(); _tmp.pop('ecolor')
axes[0].plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **_tmp)
# axes[1].errorbar(ophdata_fit.coord.l.degree, ophdata_fit.coord.distance.to(u.kpc).value,
# ophdata_fit.coord_err['distance'].to(u.kpc).value, **data_style)
axes[1].errorbar(ophdata_fan.coord.l.degree, ophdata_fan.coord.distance.to(u.kpc).value,
ophdata_fan.coord_err['distance'].to(u.kpc).value, **data_b_style)
# axes[2].errorbar(ophdata_fit.coord.l.degree, ophdata_fit.veloc['vr'].to(u.km/u.s).value,
# ophdata_fit.veloc_err['vr'].to(u.km/u.s).value, **data_style)
axes[2].errorbar(ophdata_fan.coord.l.degree, ophdata_fan.veloc['vr'].to(u.km/u.s).value,
ophdata_fan.veloc_err['vr'].to(u.km/u.s).value, **data_b_style)
# Axis limits
axes[0].set_xlim(9.,2)
axes[0].set_ylim(26.5, 33.5)
axes[1].set_ylim(5.3, 9.7)
axes[2].set_ylim(225, 325)
# Text
axes[0].set_title(name_map[name], fontsize=20)
axes[2].set_xlabel("$l$ [deg]", fontsize=18)
if i == 0:
axes[0].set_ylabel("$b$ [deg]", fontsize=18)
axes[1].set_ylabel(r"$d_\odot$ [kpc]", fontsize=18)
axes[2].set_ylabel(r"$v_r$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
axes[0].set_yticks([26,28,30,32,34])
fig.tight_layout()
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.87, 0.125, 0.02, 0.8])
fig.colorbar(cs, cax=cbar_ax)
cbar_ax.axes.set_ylabel('Stripping time [Gyr]')
fig.savefig(os.path.join(plotpath, "mockstream{}.pdf".format(k)), rasterized=True, dpi=400)
fig.savefig(os.path.join(plotpath, "mockstream{}.png".format(k)), dpi=400)
In [ ]:
for k,name_group in enumerate([all_names[:5],all_names[5:]]):
fig,allaxes = pl.subplots(2, 5, figsize=(10,5.2), sharex='col', sharey='row')
for i,name in enumerate(name_group):
axes = allaxes[:,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot, streams, prog, release_t = get_potential_stream_prog(name)
release_t = release_t/1000.
stream = streams[:,best_ix]
model_c,_model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# remove progenitor orbit
model_c = model_c[1:]
model_v = []
for v in _model_v:
model_v.append(v[1:])
# cut to just points in the window
ix = (model_c.l > 0*u.deg) & (model_c.l < 11*u.deg) & (model_c.b > 25.5*u.deg) & (model_c.b < 34.5*u.deg)
# plot model points
axes[0].scatter(model_c.l[ix], galactic.decompose(model_v[0][ix]), c=release_t[ix], **scatter_kw)
axes[1].scatter(model_c.l[ix], galactic.decompose(model_v[1][ix]), c=release_t[ix], **scatter_kw)
# Axis limits
axes[0].set_xlim(9.,2)
axes[0].set_ylim(-12,-2)
axes[1].set_ylim(-2,8)
# Text
axes[0].set_title(name_map[name], fontsize=20)
axes[1].set_xlabel("$l$ [deg]", fontsize=18)
if i == 0:
axes[0].set_ylabel(r"$\mu_l$ [${\rm mas}\,{\rm yr}^{-1}$]", fontsize=18)
axes[1].set_ylabel(r"$\mu_b$ [${\rm mas}\,{\rm yr}^{-1}$]", fontsize=18)
fig.tight_layout()
# fig.savefig(os.path.join(plotpath, "mockstream-pm{}.pdf".format(k)), rasterized=True, dpi=400)
# fig.savefig(os.path.join(plotpath, "mockstream-pm{}.png".format(k)), dpi=400)
In [ ]:
# contour stuff
levels = np.array([-2,-1,0,1]) - 0.1
fig,allaxes = pl.subplots(2, 5, figsize=(9,4.7), sharex=True, sharey=True)
for k,name_group in enumerate([all_names[:5],all_names[5:]]):
for i,name in enumerate(name_group):
ax = allaxes[k,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot, streams, prog, _ = get_potential_stream_prog(name)
stream = streams[:,best_ix]
model_c,model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# cut to just points in the window
ix = (model_c.l > 0*u.deg) & (model_c.l < 11*u.deg) & (model_c.b > 25.5*u.deg) & (model_c.b < 34.5*u.deg)
# density in sky coords
sky_ix = ix & (model_c.distance > 4.5*u.kpc) & (model_c.distance < 10*u.kpc)
grid,log_dens = surface_density(model_c[sky_ix], bandwidth=0.2)
cs = ax.contourf(grid[:,0].reshape(log_dens.shape), grid[:,1].reshape(log_dens.shape), log_dens,
levels=levels, cmap='magma_r')
# l,b
ax.plot([3.8,3.8], [31.5,32.5], **line_style)
ax.plot([5.85,5.85], [30.,31], **line_style)
# ax.set_aspect('equal')
# the data
# _tmp = data_style.copy(); _tmp.pop('ecolor')
# axes[0].plot(ophdata_fit.coord.l.degree, ophdata_fit.coord.b.degree, **_tmp)
_tmp = data_b_style.copy(); _tmp.pop('ecolor')
ax.plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **_tmp)
# Axis limits
ax.set_xlim(9.,2)
ax.set_ylim(26.5, 33.5)
# Text
ax.set_title(name_map[name], fontsize=20)
if i == 2 and k == 1:
ax.set_xlabel("$l$ [deg]", fontsize=18)
if i == 0:
ax.set_ylabel("$b$ [deg]", fontsize=18)
fig.tight_layout()
fig.savefig(os.path.join(plotpath, "mockstream-density.pdf"), rasterized=True, dpi=400)
fig.savefig(os.path.join(plotpath, "mockstream-density.png"), dpi=400)
In [ ]:
fig,allaxes = pl.subplots(2, 5, figsize=(9,5), sharex=True, sharey=True)
for k,name_group in enumerate([all_names[:5],all_names[5:]]):
for i,name in enumerate(name_group):
ax = allaxes[k,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot, streams, prog, _ = get_potential_stream_prog(name)
stream = streams[:,best_ix]
model_c,model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# ix = distance_ix(model_c.l, model_c.distance) & (model_c.l > 0) & (model_c.l < 10*u.deg) | True
ix = np.ones(model_c.l.size).astype(bool)
ax.plot(stream.pos[0][ix], stream.pos[2][ix], rasterized=True, **style)
lbd_corners = np.vstack(map(np.ravel, np.meshgrid([1.5,9.5],[26.,33.],[5.3,9.7])))
corners_g = coord.Galactic(l=lbd_corners[0]*u.degree, b=lbd_corners[1]*u.degree, distance=lbd_corners[2]*u.kpc)
corners = corners_g.transform_to(galactocentric_frame).cartesian.xyz.value
# corners = corners.reshape(3,2,2,2)[...,0].reshape(3,4)
# only pick 4 of the points to plot
corners = corners[:,corners[1] > 0.5][[0,2]]
# compute centroid
cent = np.mean(corners, axis=1)
# sort by polar angle
corner_list = corners.T.tolist()
corner_list.sort(key=lambda p: np.arctan2(p[1]-cent[1],p[0]-cent[0]))
# plot polyline
poly = mpl.patches.Polygon(corner_list, closed=True, fill=False, color='#2166AC', linestyle='-')
ax.add_patch(poly)
# ax.plot(corners[0], corners[2], ls='none', marker='o', markersize=4)
# ax.plot(-galactocentric_frame.galcen_distance, 0., marker='o', color='y', markersize=8)
ax.text(-galactocentric_frame.galcen_distance.value, 0., r"$\odot$", fontsize=18, ha='center', va='center')
# Text
ax.set_title(name_map[name], fontsize=20)
if k == 1:
ax.set_xlabel("$x$ [kpc]", fontsize=18)
# break
# break
# Axis limits
ax.set_xlim(-10,5)
ax.set_ylim(-7.5,7.5)
ax.xaxis.set_ticks([-8,-4,0,4])
ax.yaxis.set_ticks([-4,0,4])
allaxes[0,0].set_ylabel("$z$ [kpc]", fontsize=18)
allaxes[1,0].set_ylabel("$z$ [kpc]", fontsize=18)
fig.tight_layout()
fig.savefig(os.path.join(plotpath, "mockstream-xyz.pdf"), dpi=400)
# fig.savefig(os.path.join(plotpath, "mockstream-xyz.png"), dpi=400)
In [ ]:
import astropy.units as u
(300*u.km/u.s / (4*u.kpc)).to(u.mas/u.yr, equivalencies=u.dimensionless_angles())
In [ ]:
from gala.observation import distance
In [ ]:
distance(12.5 - (-2.5))
In [ ]:
def normalize_to_total_number(name, noph_stars=500): # from Bernard paper
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot, streams, prog, _ = get_potential_stream_prog(name)
stream = streams[:-1000,best_ix] # cut off at time of disruption, -250 Myr (2 stars, 0.5 Myr release)
# stream = streams[:,best_ix]
model_c,model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
ix = (model_c.distance > 4.5*u.kpc) & (model_c.distance < 15*u.kpc)
# how many points are between dense part of stream?
ix2 = ix & (model_c.l > 3.8*u.deg) & (model_c.l < 5.85*u.deg) & (model_c.b > 29.5*u.deg) & (model_c.b < 32.5*u.deg)
print(ix2.sum())
every = int(round(ix2.sum() / float(noph_stars)))
return model_c[ix]#[::every]
In [ ]:
xbounds = (2, 9)
ybounds = (27, 34)
area = (xbounds[1]-xbounds[0]) * (ybounds[1]-ybounds[0])
In [ ]:
static_gal = normalize_to_total_number("static_mw")
static_lb = np.vstack((static_gal.l.degree, static_gal.b.degree))
In [ ]:
bar_gal = normalize_to_total_number("barred_mw_8")
bar_lb = np.vstack((bar_gal.l.degree, bar_gal.b.degree))
In [ ]:
from scipy.ndimage import gaussian_filter
In [ ]:
# ddeg = (6*u.arcmin).to(u.degree).value
ddeg = (10*u.arcmin).to(u.degree).value
xbins = np.arange(xbounds[0], xbounds[1]+ddeg, ddeg)
ybins = np.arange(ybounds[0], ybounds[1]+ddeg, ddeg)
H_static,_,_ = np.histogram2d(static_lb[0], static_lb[1], bins=(xbins, ybins))
H_bar,_,_ = np.histogram2d(bar_lb[0], bar_lb[1], bins=(xbins, ybins))
bg_density = 1500 # stars / deg^2
n_bg = int(bg_density*area)
# bg_l = np.random.uniform(xbounds[0], xbounds[1], size=int(bg_density*area))
# bg_b = np.random.uniform(ybounds[0], ybounds[1], size=int(bg_density*area))
# bg_im = np.random.poisson(size=H_static.shape)
# bg_im = bg_im * float(n_bg) / bg_im.sum()
In [ ]:
bg_im = np.random.poisson(lam=42, size=H_static.shape)
bg_im.sum(), n_bg
In [ ]:
((1500 / u.degree**2) * (10*u.arcmin)**2).decompose()
In [ ]:
static_im = H_static + bg_im
bar_im = H_bar + bg_im
In [ ]:
_ = pl.hist(bg_im.ravel())
In [ ]:
# H = gaussian_filter(H, sigma=ddeg)
from scipy.stats import scoreatpercentile
_flat = static_im.ravel()
pl.hist(static_im.ravel(), bins=np.linspace(0, 1.5*scoreatpercentile(_flat,75),32), alpha=0.5);
pl.hist(bar_im.ravel(), bins=np.linspace(0, 1.5*scoreatpercentile(_flat,75),32), alpha=0.5);
In [ ]:
vmin = scoreatpercentile(_flat,2)
vmax = scoreatpercentile(_flat,98)
imshow_kw = dict(extent=[xbins.min(), xbins.max(), ybins.min(), ybins.max()],
origin='bottom', cmap='Greys', vmin=vmin, vmax=vmax, interpolation='nearest')
In [ ]:
line_style = {
"marker": None,
"linestyle": '--',
"color": 'k',
"linewidth": 1.5,
"dashes": (3,2)
}
In [ ]:
# H,_,_ = np.histogram2d(all_static[0], all_static[1], bins=(xbins, ybins))
# H = gaussian_filter(H, sigma=ddeg)
pl.imshow(static_im.T, **imshow_kw)
pl.plot([3.8,3.8], [31.5,32.5], **line_style)
pl.plot([5.85,5.85], [30.,31], **line_style)
pl.xlim(xbins.max(), xbins.min())
pl.ylim(ybins.min(), ybins.max())
pl.xlabel("$l$ [deg]")
pl.ylabel("$b$ [deg]")
pl.title("static")
# pl.tight_layout()
pl.gca().set_aspect('equal')
In [ ]:
# H,_,_ = np.histogram2d(all_bar[0], all_bar[1], bins=(xbins, ybins))
# H = gaussian_filter(H, sigma=ddeg)
pl.imshow(bar_im.T, **imshow_kw)
pl.plot([3.8,3.8], [31.5,32.5], **line_style)
pl.plot([5.85,5.85], [30.,31], **line_style)
pl.xlim(xbins.max(), xbins.min())
pl.ylim(ybins.min(), ybins.max())
pl.xlabel("$l$ [deg]")
pl.ylabel("$b$ [deg]")
pl.title("bar8")
pl.tight_layout()
In [ ]:
_density_plot_cache = dict()
In [ ]:
line_style = {
"marker": None,
"linestyle": '-',
"color": 'k',
"linewidth": 2,
}
fig,allaxes = pl.subplots(2, 5, figsize=(9,5), sharex=True, sharey=True)
vmin = vmax = None
for i,name in enumerate(all_names):
print(name)
ax = allaxes.flat[i]
if name not in _density_plot_cache:
gal = normalize_to_total_number(name)
lb = np.vstack((gal.l.degree, gal.b.degree))
H,_,_ = np.histogram2d(lb[0], lb[1], bins=(xbins, ybins))
im = H + bg_im
_density_plot_cache[name] = im.T
im = _density_plot_cache[name]
if vmin is None:
_flat = im.ravel()
vmin = scoreatpercentile(_flat,2)
vmax = scoreatpercentile(_flat,98)
imshow_kw = dict(extent=[xbins.min(), xbins.max(), ybins.min(), ybins.max()],
origin='bottom', cmap='Greys', vmin=vmin, vmax=vmax, interpolation='nearest')
ax.imshow(im, **imshow_kw)
ax.plot([3.8,3.8], [31.5,32.5], **line_style)
ax.plot([5.85,5.85], [30.,31], **line_style)
# Text
ax.set_title(name_map[name], fontsize=20)
if i > 4:
ax.set_xlabel("$l$ [deg]", fontsize=18)
ax.set(adjustable='box-forced', aspect='equal')
ax.set_xlim(xbins.max(), xbins.min())
ax.set_ylim(ybins.min(), ybins.max())
allaxes[0,0].set_ylabel("$b$ [deg]", fontsize=18)
allaxes[1,0].set_ylabel("$b$ [deg]", fontsize=18)
fig.tight_layout()
fig.savefig(os.path.join(plotpath, "densitymaps.pdf"))
fig.savefig(os.path.join(plotpath, "densitymaps.png"), dpi=400)
In [ ]:
pick = [1,8,9]
In [ ]:
style = {
"linestyle": "none",
"marker": 'o',
"markersize": 2,
"alpha": 0.3,
"color": "#555555"
}
line_style = {
"marker": None,
"linestyle": '--',
"color": 'k',
"linewidth": 1.5,
"dashes": (3,2)
}
data_style = dict(marker='o', ms=4, ls='none', ecolor='#333333', alpha=0.75, color='k')
data_b_style = data_style.copy()
data_b_style['marker'] = 's'
# contour stuff
levels = np.array([-2,-1,0,1]) - 0.1
fig,allaxes = pl.subplots(2, 4, figsize=(12,6), sharex='col', sharey='row')
for i,name in enumerate(['static_mw']+['barred_mw_'+str(j) for j in pick]):
axes = allaxes[:,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot = op.load_potential(name)
Omega = (pot.parameters['bar']['Omega']/u.Myr).to(u.km/u.s/u.kpc).value
if i == 0:
title = 'no bar'
else:
title = r'$\Omega_p={:d}\,{{\rm km/s/kpc}}$'.format(int(Omega))
pot, streams, prog, _ = get_potential_stream_prog(name)
stream = streams[:,best_ix]
model_c,model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# cut to just points in the window
ix = (model_c.l > 0*u.deg) & (model_c.l < 11*u.deg) & (model_c.b > 25.5*u.deg) & (model_c.b < 34.5*u.deg)
# density in sky coords
sky_ix = ix & (model_c.distance > 4.5*u.kpc) & (model_c.distance < 10*u.kpc)
grid,log_dens = surface_density(model_c[sky_ix], bandwidth=0.2)
cs = axes[0].contourf(grid[:,0].reshape(log_dens.shape), grid[:,1].reshape(log_dens.shape), log_dens,
levels=levels, cmap='magma_r')
axes[1].plot(model_c.l[ix], galactic.decompose(model_v[2][ix]), rasterized=True, **style)
# l,b
axes[0].plot([3.8,3.8], [31.5,32.5], **line_style)
axes[0].plot([5.85,5.85], [30.,31], **line_style)
# l,vr
axes[1].plot([3.8,3.8], [276,292], **line_style)
axes[1].plot([5.85,5.85], [284,300], **line_style)
# the data
_tmp = data_b_style.copy(); _tmp.pop('ecolor')
axes[0].plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **_tmp)
axes[1].errorbar(ophdata_fan.coord.l.degree, ophdata_fan.veloc['vr'].to(u.km/u.s).value,
ophdata_fan.veloc_err['vr'].to(u.km/u.s).value, **data_b_style)
# Axis limits
axes[0].set_xlim(9.,2)
axes[0].set_ylim(26.5, 33.5)
axes[1].set_ylim(225, 325)
# Text
axes[0].set_title(title, fontsize=18)
axes[1].set_xlabel("$l$ [deg]", fontsize=18)
if i == 0:
axes[0].set_ylabel("$b$ [deg]", fontsize=18)
axes[1].set_ylabel(r"$v_r$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
fig.tight_layout()
In [ ]:
_density_plot_cache2 = dict()
In [ ]:
line_style = {
"marker": None,
"linestyle": '-',
"color": 'k',
"linewidth": 2,
}
fig,allaxes = pl.subplots(1, 4, figsize=(11,3), sharex=True, sharey=True)
vmin = vmax = None
for i,name in enumerate(['static_mw']+['barred_mw_'+str(j) for j in pick]):
print(name)
ax = allaxes.flat[i]
pot = op.load_potential(name)
Omega = (pot.parameters['bar']['Omega']/u.Myr).to(u.km/u.s/u.kpc).value
if i == 0:
title = 'no bar'
else:
title = r'$\Omega_p={:d}\,{{\rm km/s/kpc}}$'.format(int(Omega))
if name not in _density_plot_cache2:
gal = normalize_to_total_number(name)
lb = np.vstack((gal.l.degree, gal.b.degree))
H,_,_ = np.histogram2d(lb[0], lb[1], bins=(xbins, ybins))
im = H + bg_im
_density_plot_cache2[name] = im.T
im = _density_plot_cache2[name]
if vmin is None:
_flat = im.ravel()
vmin = scoreatpercentile(_flat,2)
vmax = scoreatpercentile(_flat,98)
imshow_kw = dict(extent=[xbins.min(), xbins.max(), ybins.min(), ybins.max()],
origin='bottom', cmap='Greys', vmin=vmin, vmax=vmax, interpolation='nearest')
ax.imshow(im, **imshow_kw)
ax.plot([3.8,3.8], [31.5,32.5], **line_style)
ax.plot([5.85,5.85], [30.,31], **line_style)
# Text
ax.set_title(title, fontsize=18)
ax.set_xlabel("$l$ [deg]", fontsize=18)
ax.set(adjustable='box-forced', aspect='equal')
ax.set_xlim(xbins.max(), xbins.min())
ax.set_ylim(ybins.min(), ybins.max())
allaxes[0].set_ylabel("$b$ [deg]", fontsize=18)
fig.tight_layout()
In [ ]:
style = {
"linestyle": "none",
"marker": 'o',
"markersize": 2,
"alpha": 0.3,
"color": "#555555"
}
line_style = {
"marker": None,
"linestyle": '--',
"color": 'k',
"linewidth": 1.5,
"dashes": (3,2)
}
data_style = dict(marker='o', ms=6, ls='none', alpha=0.9, color='#2166AC')
data_b_style = data_style.copy()
data_b_style['marker'] = 's'
# contour stuff
levels = np.array([-2,-1,0,1]) - 0.1
fig,allaxes = pl.subplots(2, 4, figsize=(12,6), sharex='col', sharey='row')
for i,name in enumerate(['static_mw']+['barred_mw_'+str(j) for j in pick]):
axes = allaxes[:,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot = op.load_potential(name)
Omega = (pot.parameters['bar']['Omega']/u.Myr).to(u.km/u.s/u.kpc).value
if i == 0:
title = 'no bar'
else:
title = r'$\Omega_p={:d}\,{{\rm km/s/kpc}}$'.format(int(Omega))
pot, streams, prog, _ = get_potential_stream_prog(name)
stream = streams[:,best_ix]
model_c,model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# cut to just points in the window
ix = (model_c.l > 0*u.deg) & (model_c.l < 11*u.deg) & (model_c.b > 25.5*u.deg) & (model_c.b < 34.5*u.deg)
# density in sky coords
# sky_ix = ix & (model_c.distance > 4.5*u.kpc) & (model_c.distance < 10*u.kpc)
sky_ix = ix
grid,log_dens = surface_density(model_c[sky_ix], bandwidth=0.2)
cs = axes[0].contourf(grid[:,0].reshape(log_dens.shape), grid[:,1].reshape(log_dens.shape), log_dens,
levels=levels, cmap='magma_r')
axes[1].plot(model_c.l[ix], galactic.decompose(model_v[2][ix]), rasterized=True, **style)
# l,b
axes[0].plot([3.8,3.8], [31.5,32.5], **line_style)
axes[0].plot([5.85,5.85], [30.,31], **line_style)
# l,vr
axes[1].plot([3.8,3.8], [276,292], **line_style)
axes[1].plot([5.85,5.85], [284,300], **line_style)
# the data
axes[0].plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **data_b_style)
axes[1].plot(ophdata_fan.coord.l.degree, ophdata_fan.veloc['vr'].to(u.km/u.s).value, **data_b_style)
# Axis limits
axes[0].set_xlim(9.,2)
axes[0].set_ylim(26.5, 33.5)
axes[1].set_ylim(200, 350)
# Text
axes[0].set_title(title, fontsize=18)
axes[1].set_xlabel("$l$ [deg]", fontsize=18)
if i == 0:
axes[0].set_ylabel("$b$ [deg]", fontsize=18)
axes[1].set_ylabel(r"$v_{\rm los}$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
fig.tight_layout()
In [ ]:
style = {
"marker": '.',
"alpha": 0.75,
"cmap": custom_cmap,
"vmin": -1000, "vmax": 0.
}
line_style = {
"marker": None,
"linestyle": '--',
"color": 'k',
"linewidth": 1.5,
"dashes": (3,2)
}
data_style = dict(marker='o', ms=6, ls='none', alpha=0.9, color='#2166AC')
data_b_style = data_style.copy()
data_b_style['marker'] = 's'
# contour stuff
levels = np.array([-2,-1,0,1]) - 0.1
fig,allaxes = pl.subplots(2, 4, figsize=(12,6), sharex='col', sharey='row')
for i,name in enumerate(['static_mw']+['barred_mw_'+str(j) for j in pick]):
axes = allaxes[:,i]
# path to file to cache the likelihoods
cache_file = os.path.join(RESULTSPATH, name, "mockstream", "ln_likelihoods.npy")
lls = np.load(cache_file)
best_ix = lls.sum(axis=1).argmax()
pot = op.load_potential(name)
Omega = (pot.parameters['bar']['Omega']/u.Myr).to(u.km/u.s/u.kpc).value
if i == 0:
title = 'no bar'
else:
title = r'$\Omega_p={:d}\,{{\rm km/s/kpc}}$'.format(int(Omega))
pot, streams, prog, release_t = get_potential_stream_prog(name)
stream = streams[:,best_ix]
model_c,_model_v = stream.to_frame(coord.Galactic, galactocentric_frame=galactocentric_frame, vcirc=vcirc, vlsr=vlsr)
# remove prog orbit
model_c = model_c[1:]
model_v = []
for v in _model_v:
model_v.append(v[1:])
# cut to just points in the window
ix = (model_c.l > -2*u.deg) & (model_c.l < 14*u.deg) & (model_c.b > 23.5*u.deg) & (model_c.b < 36.5*u.deg)
axes[0].scatter(model_c.l[ix], model_c.b[ix], rasterized=True, c=release_t[ix], **style)
axes[1].scatter(model_c.l[ix], galactic.decompose(model_v[2][ix]), rasterized=True,
c=release_t[ix], **style)
# l,b
axes[0].plot([3.8,3.8], [31.5,32.5], **line_style)
axes[0].plot([5.85,5.85], [30.,31], **line_style)
# l,vr
axes[1].plot([3.8,3.8], [276,292], **line_style)
axes[1].plot([5.85,5.85], [284,300], **line_style)
# the data
axes[0].plot(ophdata_fan.coord.l.degree, ophdata_fan.coord.b.degree, **data_b_style)
axes[1].plot(ophdata_fan.coord.l.degree, ophdata_fan.veloc['vr'].to(u.km/u.s).value, **data_b_style)
# Axis limits
axes[0].set_xlim(9.,2)
axes[0].set_ylim(26.5, 33.5)
axes[1].set_ylim(200, 350)
# Text
axes[0].set_title(title, fontsize=18)
axes[1].set_xlabel("$l$ [deg]", fontsize=18)
if i == 0:
axes[0].set_ylabel("$b$ [deg]", fontsize=18)
axes[1].set_ylabel(r"$v_{\rm los}$ [${\rm km}\,{\rm s}^{-1}$]", fontsize=18)
fig.tight_layout()
In [ ]:
from ophiuchus.experiments import LyapunovGrid
In [ ]:
bins = np.linspace(0.3,1.1,12)*u.Gyr
fig,axes = pl.subplots(1, 4, figsize=(12,3.5), sharex='col', sharey='row')
for i,name in enumerate(['static_mw']+['barred_mw_'+str(j) for j in pick]):
pot = op.load_potential(name)
Omega = (pot.parameters['bar']['Omega']/u.Myr).to(u.km/u.s/u.kpc).value
if i == 0:
title = 'no bar'
else:
title = r'$\Omega_p={:d}\,{{\rm km/s/kpc}}$'.format(int(Omega))
gr = LyapunovGrid.from_config(os.path.join(RESULTSPATH,name,"lyapunov"),
os.path.join(RESULTSPATH,"global_lyapunov.cfg"),
potential_name=name)
d = gr.read_cache()
ftmle = (d['mle_avg']*1/u.Myr)
lyap_time = (1/ftmle).to(u.Myr)
axes[i].hist(lyap_time, bins=bins.to(u.Myr))
axes[i].set_xlim(200,1300)
axes[i].set_ylim(0,60)
axes[i].xaxis.set_ticks([300,600,900,1200])
axes[i].yaxis.set_ticks([0,20,40,60])
axes[i].set_xlabel(r"$t_\lambda$ [Myr]", fontsize=18)
axes[i].set_title(title, fontsize=18)
axes[0].set_ylabel('$N$')
fig.tight_layout()